import numpy as np
import matplotlib.pyplot as plt
from scipy.special import factorial

plt.rcParams['font.sans-serif'] = ['SimHei']  # 用黑体显示中文
plt.rcParams['axes.unicode_minus'] = False    # 正常显示负号

# 参数化圆周 |z| = 2
def circle_path(t, radius=2):
    return radius * np.exp(1j * t)

# 圆周 |z| = 2 的参数方程的导数
def derivative_circle_path(t, radius=2):
    return 1j * radius * np.exp(1j * t)

# 被积函数 f(z) = e^z / z^n
def integrand(z, n):
    return np.exp(z) / (z**n)

# 计算积分
def cauchy_integral(n, radius=2, n_points=100):
    t = np.linspace(0, 2*np.pi, n_points) # 参数区间上的离散值
    z = circle_path(t, radius) # 积分路径上的离散值
    dz = derivative_circle_path(t, radius) * (t[1] - t[0]) # 积分路径上的微分
    
    integral = np.sum(integrand(z, n) * dz) # 黎曼和
    return integral

# Cauchy积分公式（高阶导数公式）: 
# 对于函数f(z)=e^z, f^(n-1)(0) = (n-1)! / (2*pi*i) * integral
# 即 integral = (2*pi*i) * f^(n-1)(0) / (n-1)!
# 因为 e^z 的各阶导数都是 e^z, 所以 e^z 在0处的各阶导数都是1
def analytical_result(n):
    if n <= 0:
        return 0  # 如果n不是正整数
    return (2 * np.pi * 1j) / factorial(n-1) if n > 0 else 0

print("验证Cauchy积分公式（高阶导数公式）:")
print("n值\t数值积分结果\t解析解\t\t误差")

for n in range(1, 6):
    numerical = cauchy_integral(n)
    analytical = analytical_result(n)
    error = abs(numerical - analytical)
    print(f"{n}\t{numerical:.4f}\t{analytical:.4f}\t{error:.2e}")

